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Abstract 

The use of adaptive mesh refinement (AMR) techniques is crucial for accurate and 
efficient simulation of higher dimensional spacetimes. In this work we develop an 
adaptive algorithm tailored to the integration of finite difference discretizations of 
wave-like equations using characteristic coordinates. We demonstrate the algorithm 
by constructing a code implementing the Einstein-Klein-Gordon system of equations 
in spherical symmetry. We discuss how the algorithm can trivially be generalized to 
higher dimensional systems, and suggest a method that can be used to parallelize 
a characteristic code. 



1 Introduction 



The investigation of the characteristic structure of a system of partial dif- 
ferential equations (PDEs) gives valuable insight on the behavior of allowed 
solutions. Analytical studies have long benefited from understanding and em- 
ploying this knowledge. In the numerical realm, efficient algorithms have been 
developed which exploit the underlying characteristic structure. These algo- 
rithms are obtained by formally transforming to characteristic variables, re- 
alizing how these need to be updated, and then employing this information 
to construct an algorithm using the original variables. However, the explicit 
integration along characteristics, where coordinates are chosen adapted to the 
characteristic directions, has received considerably less attention. Indeed, this 
option has so far only been actively pursued within general relativity (GR), 
although as discussed in [1,2] there are powerful reasons to consider this option 
in systems described by wave-like equations. 
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The characteristic formulation of general relativity has since its introduction 
in the 60's played an important role as a tool to investigate different aspects of 
the theory (see for instance [3,4,5,6,7,8].) A clean picture of the effect of grav- 
itational waves and their geometric manifestation, links between the structure 
of future null infinity and some interior sources, and the analysis of singularity 
structure are some areas that have effectively been tackled by this approach. 

In recent years, the formalism has also displayed its usefulness in the numerical 
arena. Investigations of critical phenomena[9,10,ll,12,13], wave propagation 
on non-trivial backgrounds[14] and simulations of spacetimes containing black 
holes[15,16,17,18] or neutron stars[19] have benefited from exploiting this ap- 
proach. Numerical simulations of black hole spacetimes are of considerable 
importance for the detection and analysis of gravitational waves that could 
be measured by the new generation of gravitational wave detectors. Among 
the systems expected to produce gravitational waves of sufficient intensity are 
those containing black holes and neutron stars, and it would be ideal if the suc- 
cess obtained simulating single black holes with 3D characteristic codes could 
be translated to these systems. Preliminary indications that this is likely the 
case for a subset of possible black hole-neutron star binaries can be found in 
[20,19], where higher dimensional, non- vacuum scenarios (yet simpler than the 
3D binary black hole-neutron star system) have been accurately simulated. 

Unfortunately, the computational requirements for modeling a binary black 
hole-neutron star system are quite high if one of the goals is to predict the 
waveforms produced by it in a quantitative manner. Since the expected energy 
output in gravitational waves is at most a few percent of the total mass of the 
system, the numerical simulation must guarantee that any systematic (numer- 
ical) error is well below this target. The major issue in achieving an accurate 
description of the system, when a stable discretization has been obtained, 
is to adequately resolve the different length scales involved. These scales are 
naturally defined by the stellar dynamics, the black hole, and the distant 
weak-field regime where the gravitational waves are extracted. Covering all of 
these scales with a uniform grid (in a finite-difference based numerical simu- 
lation) of sufficient resolution to accurately model the smallest features is not 
only a waste of resources, but may be impossible to achieve on contemporary 
computers. Therefore, we need to use techniques that better exploit available 
resources. One such technique is adaptive mesh refinement (AMR), whereby 
the simulation can dynamically adjust the grid resolution in different regions 
of the domain to adequately resolve all features with sufficient, but not excess, 
resolution. 

In this paper we investigate the use of AMR techniques for characteristic evo- 
lution. Although we concentrate on the GR case, the algorithms presented here 
can readily be applied to many other equations. The use of characteristic AMR 
in GR has been partially addressed by several authors in the past[9,10,12]. 
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however, in those cases the algorithms were geared to studying gravitational 
collapse and singularity structure in spherical symmetry and the techniques 
do not generalize to higher dimensional systems. Here we present an AMR 
algorithm for characteristic evolution that can be applied to scenarios with 
an arbitrary number of non-trivial spatial dimensions. Furthermore, the algo- 
rithm does not place significant restrictions on the discretization scheme, and 
hence existing unigrid codes (as described, for instance, in [21,15,22,23,17,18]) 
can, in principle, be incorporated into the adaptive framework in a straight- 
forward manner. To demonstrate the basic algorithm, and show that it can 
adapt to dynamical features of a solution, we have implemented a spherically 
symmetric code to solve the Einstein-Klein-Gordon system. 

The rest of the paper is organized as follows. In Sec. 2 we describe the coor- 
dinate system, set of equations and corresponding discretization scheme that 
we will employ in the example problem. We use a coordinate system with a 
single null direction, however, the AMR algorithm is most easily presented 
for double null coordinates; hence in Sec. 3 we describe the adaptive scheme 
for double null evolution. In section 4 we mention how the basic algorithm is 
extended to a coordinate system with a single null coordinate, how additional 
non-trivial spacelike dimensions can be handled, mention some problem spe- 
cific features such as the numerical dissipation and interpolation operators we 
use, and describe a method that can be used to parallelize a characteristic 
code (with or without AMR). In Sec. 5 we present results from our spherically 
symmetric code, and give some concluding remarks in Sec. 6. 



2 The spherically symmetric Einstein-Klein- Gordon problem in 
the chciracteristic approach 

Einstein's equations can be expressed in notational form as Gab — '^'^Tab, 
where Gab is the Einstein tensor and Tab the stress energy tensor of the matter 
distribution. In the particular case of a scalar field $ coupled to GR, Tab results 



where gab is the metric tensor, and m is the mass parameter of the scalar field. 

We introduce a coordinate system adapted to incoming null hypersurfaces in 
the following way: incoming lightlike hypersurfaces are labeled with a param- 
eter V, each null ray on a specific hypersurface is labeled with [9, 0) and r is 
introduced as a surface area coordinate (i.e. surfaces at r = const, have area 
47rr^). In these coordinates, the metric takes the Bondi-Sachs form [3, 4] 



in[24] 



V„$V6$ - 1/2 Qab (Vc$V^$ + m'$') 
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(2) 



Note that the choice of incoming null surfaces is merely for convenience; the 
trivial change /3 — >■ /9 + i7r takes the line element (2) into the one corresponding 
to outgoing null surfaces (further details can be found in [25]). The algorithm 
presented here can easily be modified to handle the outgoing case. 

In order to express the equations of motion in a simpler form, we introduce 
the following variables 



V = -r + g (3) 
* = r$ . (4) 

The resulting equations (provided by R^r = 87r$^^, Rqq = Stt^'^q and □$ = 
respectively, where Rab is the Ricci tensor) reduce to: 



A.=2.^^:%^, (5) 



g^r = l- e^^ (l - 47rm^*2) , (6) 

m,r. = - (i - f ) + ) ^ (*,^ - 7) + e''^'* (7) 

A properly posed problem requires data to be given on an initial v = hypcr- 
surface (which consists only of the unconstrained ^') and consistent boundary 
data on an intersecting surface (which includes ^, g and /3). The boundary 
data comprises gauge, physical and constrained data. For instance, given the 
value of /? and \E' on an r = const, surface the value of g is determined by 
the remaining Einstein equation R^^ = 87i^\ (which we call the consistency 
equation). Next we describe the particular setting used in the tests performed 
throughout this work. 



2.1 Coordinate conditions and boundary data 



For our present purposes we choose the outer boundary to coincide with past 
null infinity (/~). At /~ we fix a Bondi coordinate system, therefore P — 
(i.e. V represents the affine time of observers at /~). Since past infinity is a 
null hypersurface, data for ^ is unconstrained and in fact is intimately related 
(just a time derivative difference) to the "Bondi News" which is the incoming 
radiation from past null infinity. We provide data for ^' arbitrarily, and then 
determine the value of g using the consistency equation, which on I~ reduces 
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to: 



5,. = 87r(vl/,,)2. (8) 
2.2 Numerical Implementation 

We discretize equations (5,6,7,8) using second order finite difference (FD) 
tecfiniques. In order to include /~ in our computational grid, we introduce a 
compactified radial coordinate x = r/{l+r) (hence r G [0, oo)) [23] . This 
coordinate is uniform in the domain [0, 1] and is discretized by Xi = {i — 1) Ax, 
with i — l,2,...,Nx; similarly, v is discretized so that v"' — {n — l)Av, with 
n — 1,2,...,N„. As is customary, we label a grid function / at coordinate 
location {xi. y") by /" — sec Fig. 1 below for a schematic representation of the 
coordinate system and discretization. The boundary data is specified along 
the pair of intersecting null surfaces x = 1 (/^) and w = 0. Integration of the 
evolution and hypersurface equations (5-7) then proceeds via a sequence of 
radial integrations (inwards, from x — 1 — Ax to x — along each of the 
null surfaces from v — Ai> to v — Vi. Specifically, the discretized hypersurface 
equations read: 



/?^' = /3^+Y-2 7rAa;^^^ Ml - x,)^ - ^if , (9) 

= - Ax (1 - x^f (1 - e^^\°{l - AT^m'm)) , (10) 

where 

Xc = {xi + Xi+i)/2, (11) 

*|, = (vl/-+i + ^-+/)/2, (12) 

*,x|. = (vl>^+/-vl>-+i)/Ax; (13) 



and similarly for l3\o- The evolution equation for ^ is evaluated at the point 
(f "^+Af /2, Xi+Aa;/2) by including the points (u""*"-^, Xj), (v""*"-^, Xj+i), (v'^"'"-^, Xi+2), 
(i;",Xi_i),(w",Xj) and (w",Xi+i): 

*r^=-(-*ri^(-c+^(i-xc)^Fi) 

+ ^1:2 (^(1 - ^cfF,) + Xe(^r+l - 

-FiAxA^;^,^|e(l-Xe) 
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,x\c 



,x\c 




(14) 



where 



9,x\c 




)/2/Ax 



(15) 
(16) 
(17) 



and analogously for P\c and "if^xlc- For each integration step (14) is first 
used to solve for then (9) is solved for /Sf"*"^, and finally wc obtain g]^~^^ 

from (10). Note that at the first radial point in from /" (at a; = 1 — Ax) there 
are insufficient points available to evaluate g^x\c via the above scheme, and so 
there we simply propagate ^ inwards using 



Regularity conditions at a; = are explicitly enforced for ^ {^(x = 0, ?;) = 0} 
and (3 {l3x{x = 0,v) — 0}. Furthermore we set ^(x = Ax,v) using fourth 
order interpolation: 



On the outgoing (,x = 1) and ingoing {v = 0) initial characteristics, we freely 
specify We set g{x = l,v = 0) = 2mo, where mo is the initial, asymptotic 
mass of the spacetime. We then integrate (9-10) to obtain P{x,v = 0) and 
g{x,v — 0), and integrate (8) along x — 1 (using a similar discretization to 
(10)) to obtain g{x — l,v). 

We employ black hole excision techniques to solve for spacetimes containing 
a black hole, and consequently a geometric singularity. Here the underlying 
assumption is that cosmic censorship holds: any singularity will be hidden 
by an event horizon (black hole), and hence cannot causally infiuence the 
region outside the black hole[26]. This feature is exploited by placing an inner 
boundary inside the event horizon, preventing the simulation from getting 
too close to the singularity. One cannot determine the location of the event 



(18) 




(19) 
(20) 
(21) 
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{x,v) = {l,0) 



Fig. 1. A schematic representation of the {x, v) coordinate system (2) and discretiza- 
tion scheme on a spactime diagram, where null directions are at 45 degree angles 
relative to the vertical, and timelike (spacelike) curves have tangent vectors less 
than (greater than) 45 degrees. The coordinate lines v = const, are ingoing null 
curves, and x = const. (< 1) are timclike curves of constant areal radius. For refer- 
ence, outgoing null curves (thinner dotted lines) are also shown on the plot. Note 
that X = const, becomes null in the limit x ^ 1, corresponding to past null infinity. 
The discretization of a variable / is also shown on the figure — the points labelled 
correspond to those that need to be provided (in general) as "initial data" to solve 
for the unknown /""''^ at the interior point {xi,v"'^^). 

horizon prior to solving for the geometry of the entire spacetime; however, 
we can use the location of an apparent horizon to tell us where to excise, 
for under reasonable assumptions the apparent horizon always lies inside the 
event horizon [2 7]. An apparent horizon is defined as the outermost, closed 
surface whose outgoing null rays form a non-divergent front (i.e the surface 
is 'trapped'). The location of the apparent horizon is given by y = in 
(2). In practice, when we detect an apparent horizon, we excise the portion 
of the grid interior to it, though we leave a small buffer zone between the 
apparent horizon and actual excision surface (which is inside the apparent 
horizon). We implement excision by extrapolating all variables, using fourth 
order extrapolation, to all points interior to the surface of excision. Thus the 
"solution" we obtain inside the black hole is not physically meaningful, but 
the excision does not adversely affect the exterior part of the solution that we 
are interested in. 
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3 AMR in double null coordinates 

Here we motivate and describe our AMR algorithm for characteristic codes. 
The sahent features of the algorithm are best demonstrated in a double null 
coordinate system; therefore we first describe the algorithm in detail for this 
case, and then discuss the modifications for the spacelike-null situation in the 
following section. 

Our AMR algorithm is modeled after the Berger and Ohger (B&O) algorithm 

[28] for hyperbolic, Cauchy problems. The B&O algorithm has several desir- 
able features that we have used as cornerstones in building the scheme for 
characteristic codes: 

• the computational domain is decomposed into a grid hierarchy, whereby 
the partial differential equations are discretized using identical unigrid finite 
difference schemes on each grid within the hierarchy. 

• dynamical regridding is performed via local truncation error (TE) estimates. 

• the recursive evolution algorithm makes efficient use of resources in both 
space and time, for the grids are always evolved with a time step set by the 
local spatial discretization scale (to satisfy the CFL condition), and not by 
the smallest scale within the problem. 

It is not possible to directly apply the B&O algorithm in a double null co- 
ordinate system by (for instance) treating one of the null coordinates as the 
"spatial" coordinate and the other as "time" , and then integrating the "spa- 
tial" surfaces in "time". This is because propagation along the null surface 
masquerading as a spatial surface will be instantaneous, and hence the local 
TE estimation scheme will not be able to track corresponding features of the 
solution. The key to adapting B&O to characteristic codes is to effectively con- 
sider each null direction as "time", and then to simultaneously evolve along 
both. 

In the double null evolution algorithm the structure of the grid hierarchy goes 
hand-in-hand with the evolution scheme, so we first describe the hierarchy in 
detail before presenting the evolution scheme. 



3.1 AMR grid hierarchy 

For the following discussion we will consider the discretization of a double null 
coordinate system {u,v), where u = const, labels an outgoing null curve, and 
V = const, an ingoing null curve. For example, a coordinate transformation 
of the form du — — e^" {Vdv/r + 2dr), with e^"^^''') an integrating factor, will 
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bring (2) into the form 

ds^ = -e^^dudv + r'^dQ^, 



(22) 



where ^ and r are now considered functions of u and v. 

The AMR grid hierarchy (sec Fig. 2 below for an example) consists of a se- 
quence of N levels: £i,£2, ■■■An-, where all grids at level m are discretized on 
a uniform mesh with [Awm, Au^] = [Awi/p"^"^, Avi/p™"^], and p is the re- 
finement ratio. Therefore, higher levels (larger level numbers) are composed of 
finer grids. In this paper we only consider p — 2, however the generalization to 
arbitrary integer values of p > 2 is straight-forward (as is the generalization 
to different refinement ratios for each level). All grids at level m are entirely 
contained within grids at level m — 1. What we mean by this is that the co- 
ordinate region spanned by the union of grids at level m is a subset of the 
coordinate region of the union of grids at level m — 1. Furthermore, we require 
that the meshes of all levels be aligned such that any point {u^,v^) in a grid 
at level m — 1 within the region of overlap between levels m and m — 1 (a 
parent point) be coincident with a point on a grid at level m (a child point). 
This alignment of grids between levels is essential for the recursive evolution 
algorithm (and useful for truncation error estimation), as will be explained 
in the next section. As an aside, note that the characteristic AMR algorithm 
could still allow for rotation of child grids in higher dimensional (> 2) sim- 
ulations in the same sense as the original B&O algorithm, however, here the 
rotation would be within a subspace orthogonal to u — const., f = const.. In 
other words, we only require rigid alignment of the two "time" coordinates u 
and V. 

3.2 Evolution Scheme 

Evolution of PDEs discretized on the kind of grid hierarchy just described 
proceeds by recursively evolving a particular sequence of unigrid unit cells, 
from the coarsest to finest levels, and then propagating the solution obtained 
on the finer levels back to the coarser ones. The unit cell for a typical double 
null discretization scheme is shown in Fig. 3. For a concrete example, consider 
the spherically symmetric wave equation on a fiat background (whose line 
element is ds"^ — —dudv + r^dQ^, where r — {v — u)/2): 



A second order accurate finite difference version of (23), discretized on the 
unit cell of Fig. 3, is[29] 



(23) 
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Composite Grid Structure 




Level 1 / 




Level 2 



Level 3 



Level 4 




Fig. 2. An example of a double null grid hierarchy. The upper-most figure shows 
the entire hierarchy, composed of four levels in this case. Each level is stored inde- 
pendently of the others; in the lower plots we show the structure of each of these 
levels. 
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Fig. 3. A fundamental unit cell of the double null evolution scheme. The unit cell con- 
sists of the four points A, B, C and D. In the plot we also schematically show part of 
the data structure we use to represent the hierarchy, namely, a set of point-structures 
(or points for short), linked together to form the mesh. Each point contains north 
(ra), south (s), east (e) and west (w) links to adjacent points at the same level, as 
shown for point C. In addition, certain points will have parent and/or child links to 
corresponding points in the parent and/or child levels (not shown in the figure). 

AuAv 

2r \ 2Au 2Av 

Initial data for must be specified on the initial ingoing and outgoing char- 
acteristics, which amounts to specifying at points B, C and D. Evolution 
to point A then proceeds by solving (24) for (pA- Evolution of (p over an en- 
tire uniform mesh of cells is then trivial (ignoring boundary conditions) — the 
unit cell evolution scheme is repeatedly applied, in arbitrary order, to all cells 
where past values (corresponding to points A, C and D) of (p are known, until 
all points in the grid are solved for. 

The extension of this unigrid evolution scheme to an adaptive hierarchy is 
for the most part straightforward, and follows the B&O scheme rather closely 
in spirit. The evolution algorithm consists of a particular sequence of single, 
unit cell evolution steps. The order in which cells arc traversed over the hi- 
erarchy is dictated by causality, in that we can only evolve to point A if the 
"initial/boundary data" at points B, C and D are known. However, notice in 
Fig. 2 (and below in the first frame of Fig. 5, depicting a sample evolution) 
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that this initial data will not be available on a finer level m > 1 with initial 
surfaces u = U„i oi v = Vm that arc interior to the computational domain 
boundaries, i.e. where 1/^ > Ui or > Vi (we assume that at u = Ui and 
V — Vi the entire initial grid structure is supphed — see Sec. 4 for a discussion 
on how the initial hierarchy could be calculated) . The solution to this problem 
is to always evolve coarser, parent cells first, and use the solution obtained on 
a parent cell to set initial data, via interpolation, at child cell points bounding 
a newly refined region. Then, evolution on the set of child cells is performed, 
recursively evolving additional finer levels if present. The solution obtained 
at points coincident with points on the parent level are injected back to the 
parent level, to maintain a single-valued solution where the most accurately 
known values are stored at all points in the hierarchy. 

In theory, using interpolation to set interior initial data of fine regions will not 
adversely affect the consistency of the FD approximation to the PDEs if the 
hierarchy is generated via local truncation error (TE) estimates ^ . For then, 
prior to the surface u = Ue {oi v = Vi) when a new fine level i is added, the 
local TE on level £ — 1 will have the same order of magnitude as the maximum 
allowed TE, and hence if a sufficiently high order of interpolation is used to 
initiahze the fields ai, u — {v — Vi) the solution error on level i will differ 
from that on level £ — 1 by an amount less that the local TE there. In practice, 
interpolation often introduces high frequency solution components ("noise") 
that produce a significant amount of error when propagated away from the 
refinement boundaries during subsequent evolution. Adding numerical dissi- 
pation to the FD scheme can significantly reduce this noise, as can the choice 
of interpolation operator and the frequency of regridding. These issues will be 
discussed in more detail later on in this section, and in Sec. 4.1.1. 

The rule that we always evolve a parent cell before any child cells is naturally 
implemented via a recursive subroutine, which is summarized in pseudo-code 
in Fig. 4 below. The steps taken in a sample, 3-lcvel evolution is depicted in 
Fig. 5. One of the differences of the characteristic AMR algorithm, compared 
to the B&O algorithm, is that a slightly more complicated data structure 
is needed to efficiently represent the dynamical hierarchy. In B&O, the grid 
hierarchy is calculated over the entire spatial hypersurface at a given time, 
and hence the grids at a given level can efficiently be stored as a list of one 
dimensional arrays (in a 1+1 D simulation). In the characteristic algorithm, 
the structure of the hierarchy is revealed point-by-point, simultaneously in the 
u and V directions as one evolves, and hence one cannot pre-allocate similar 
one dimensional arrays. We have chosen to use a data structure where a point 



We assume that the solution and truncation error estimates are sufficiently 
smooth functions of u and v. However, in principle one can specify non-smooth, 
though continuous initial data in (p on v = vq and u = uq, together with the appro- 
priate initial hierarchy. 
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(m*, ?,'•'), at some level £, is the fundamental unit of data. The mesh is then con- 
structed by linking together adjacent points at the same level with north(n), 
south(s), east(e) and west(w) pointers as depicted in Fig. 3, and linking points 
at the same coordinate location in levels i — 1 and i + 1 via parent and child 
pointers respectively. In the pseudo-code in Fig. 4 we have used C program- 
ming language notation to represent links; for example, referring to Fig. 3, 
A = B ^ n, and B = A — ^ s. In the following paragraphs we will discuss key 
lines of the pseudo-code in detail. 

The function evolve-uniLcell(C) listed in Fig. 4 takes, if possible, a single 
evolution step to the causal future of the point C, at the level of point C. 
Prior to solving the PDEs in line 16, the hierarchy is extended to points A, B 
and D if necessary (lines 6 — 14). As discussed in the preceding paragraphs, 
when the program execution reaches line 6, points B and D will always exist 
on the coarsest level of the hierarchy, and at interior points of all levels; only 
on the boundaries of a refined region might one need to create these points 
and initialize them via interpolation from the parent grid. Once the equations 
have been solved at point A, if the TE at C is greater than the maximum 
allowed, we recursively evolve the four unit cells of the child level that occupy 
the same region as the unit cell of point C (fines 18 — 23) ^ . Note that because 
of the manner in which we refine (lines 31 — 34), C will always have a child 
point in the hierarchy if the TE at C is greater than the maximum allowed 
value. After all the child-levels have been evolved to point A, the solution 
obtained there on the finest level is injected back to the current level (line 30) ; 
i.e. we replace the variables at A with those at the child of A (the recursion 
guarantees that the values stored at the child of A will have come from the 
solution obtained on the finest level of the hierarchy containing A) . 

We compute the truncation error estimate in lines 25 — 27 using a self-shadow 
hierarchy technique [30], which is a variant of a shadow hierarchy. In the 
traditional B&O algorithm, when truncation error estimates are needed at 
regridding time two copies of the subset of the hierarchy over-which regrid- 
ding will be performed is made; the first is an identical copy of the hierarchy, 
while the second is a 2 : 1 coarsened version of the first. Then a single evo- 
lution step is taken on the coarsened copy, and two evolution steps (with the 
same Courant factor) arc taken on the fine copy, i.e. the copies are evolved 
to the same coordinate time. The TE is then computed as a point-wise norm 
of the difi'crence between the solutions obtained on the two copies (which are 
subsequently deleted). A shadow-hierarchy economizes this process by always 
evolving a 2:1 coarsened version (the "shadow") of the main hierarchy in con- 
junction with the main hierarchy. The solution obtained on the main hierarchy 
is then periodically injected into corresponding grids of the shadow hierarchy. 



^ In general, a refinement ratio of n : 1 will require a sequence of evolution steps 
on the child level at this stage in the algorithm. 



13 



1: logical fionction evolve_unit_cell (point C) 

2: if (can_take_step(C)==f alse) then 

3: return(f alse) ; 

4: end if 

5: 

6: if (point B=C->w does not exist) then 

7: create B and link it into the surrounding mesh; 

8: initialize variables at B via interpolation from parent points; 

9: end if 

10: if (point D=C->n does not exist) then 
11: create D and link it into the surrounding mesh; 

12: initialize variables at D via interpolation from parent points; 

13: end if 

14: create point A=D->w(=B->n) and link it into the surrounding mesh; 
15: 

16: solve for the variables at A via the evolution equations; 
17: 

18: if (TE(C)>maximum_TE) then 

19: evolve_unit_cell(C->child) ; 

20: evolve_unit_cell(C->child->w) ; 

21: evolve_unit_cell(C->child->n) ; 

22: evolve_unit_cell(C->child->n->w) ; 

23: end if 

24: 

25: if (A->parent exists) then 
26: compute_TE(A) ; 

27: end if 
28: 

29: if (A->child exists) then 
30: inject variables from A->child to A; 

31: else if (TE(A) > maximum_TE) then 
32: create point A->child; 

33: initialize variables of A->child with the corresponding values at A; 

34: end if 

35: 

36: (one can safely delete points to the causal past of C here) 
37: 

38: return (true ) ; 

39: end of function evolve_unit_cell; 

40: 

41: function evolve_hierarchy () 

42: Nul=number of points in the u-direction in the base level(l); 

43: Nvl=number of points in the v-direction in the base level; 
44: 

45: create the hierarchy corresponding to the initial surfaces 
46: u=uO and v=vO, and initialize the variables there; 

47: 

48: set point Ci = grid point at (uO.vO) on the base level; 

49: for i=l to (Nvl-1) do 

50: Cj=Ci; 

51: for j=l to (Nul-1) do 

52: evolve_unit_cell(Cj) ; 

53: Cj=Cj->w; 

54: end do 

55: Ci=Ci->n; 

56 : end do 

57: end of function evolve_hierarchy ; 



Fig. 4. A pseudo-code description of the adaptive evolution algorithm. The function 
evolve-unit-cell(C) recursively evolves the PDEs on all points in the grid hierarchy 
at the level of point C and higher one unit cell (of the level of C) to the causal 
future of C (i.e. to point A in Fig. 3). The function evolve Jiierarchy() demonstrates 
one possible sequence of evolvejunit_cellQj calls that can be used to solve the PDEs 
over the entire hierarchy (we arbitrarily decided to evolve in u first, then v). 
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Fig. 5. Steps in a sample, 3-level evolution. Frame 1 shows the initial hierarchy. 
The subsequent frames show how the hierarchy is recursively created in a single 
coarse grid (at level 1) evolution step. An arrow in a frame indicates which grid 
point is being updated during that step. Frames 6 and 7 show the only steps in 
this example where data needs to be interpolated from a parent level (at the point 
to the "south" of the arrowed-point in each case). After step 9, we assume that 
level 3 is unrefined, and hence the final evolution step depicted in frame 10 takes 
place at level 2. Note that in the algorithm described in the paper, unrefinement 
(and equivalently refinement in frame 6) would not occur here; using a self-shadow 
hierarchy, refinement/unrefinement of level 3 can only occur at points where level 
2 and 1 are in sync. For brevity we ignore this aspect of the algorithm here. 
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A self-shadow hierarchy further economizes the truncation error estimation 
process by noting that, due to the recursive nature of the evolution algorithm, 
information required to compute a TE at level i is "naturally" available prior 
to the injection step from level i — 1 to £. For at that stage in the evolu- 
tion process, the solution obtained at the common point A on both levels has 
been calculated via independent evolution at two discretization scales, start- 
ing from identical "initial data" one unit cell to the past of point A on level 
i — 1. Thus, metaphorically speaking, a hierarchy can act as its own shadow 
at points where there are at least two levels of refinement. To implement a 
self-shadow hierarchy requires that the base level (level 1) always be fully re- 
fined, and we then define the TE at a point in level £, £ > 1, via the difference 
in the solutions obtained there and at the corresponding point in level £ — 1, 
prior to injection from level £ to £ — 1. 



3.2.1 More on Truncation Error Estimates 

We conclude this section by discussing a few practical details concerning the 
computation of the TE. The point-wise TE computed using solutions to wave- 
like finite-difference equations is in general oscillatory in nature, and will tend 
to go to zero at certain points within the computational domain (we will dis- 
cuss this in more detail in the next paragraph), even in regions of relatively 
high truncation error. We do not want such isolated points of (anomalously) 
small TE to cause temporary unrefinement, for experience suggests that re- 
finement boundaries are often a significant source of unwanted high-frequency 
solution components. Even though we can, to some degree, eliminate the high- 
frequency components via dissipation (see Sec. 4) one would like to avoid sit- 
uations that produce "noise" as much as possible. Therefore, in practice, the 
TE we use to determine whether we refine or unrefine at a given point is an 
average of the point- wise TE taken over several cells to the past of the point. 
Also, note that when using a self-shadow hierarchy, we only compute a point- 
wise TE when point A is in sync with its parent; we then define the TE at 
the three points A ^ n, A ^ w and ^ ^ n — > w to be identical to that of A. 

To give a more quantitative description of the nature of the TE, we will analyze 
the sample wave equation and discretization given above (23-24). Decompose 
a solution to the finite difference equation C4> = (24) as 

(/) = 00 + 0e, (25) 

and decompose the difference operator C as 

C^Co + JCe, (26) 



where jCq — d^v + (l/2r)(9u — d^) is the continuum operator (23), and 0o 
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satisfies the continuum wave equation {Co<Po — 0). Hence 0e is the truncation 
error. For the discretization in (24), the operator Cg takes the form 
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where h denotes either Au or Av. Then 

= (£o + >Ce)(0O + 0e) 
£o0e + >Ce0O, 



(27) 



(28) 



where in the last step we have assumed that the truncation error 0e is of order 
h^, and so to leading order we can ignore the term >Ce0e- Considering equation 
(28) to be an evolution equation for the truncation error 0e, ^-nd assuming 
that 00 is given, we can see that 0e satisfies the continuum wave equation 
with source term —Ce(t>Q'- 



>Co0e ~ ->Ce0o- 



(29) 



Therefore, from (27), the leading order part of the truncation error will be 
proportional to third and fourth derivatives of 0o. During a numerical evo- 
lution, if we are in the convergent regime, then a truncation error estimate 
computed as described in the preceding paragraphs will be a good approxi- 
mation to the actual truncation error 0e, and the numerical will be close to 
the desired continuum solution 0o; hence we would expect the TE estimate 
to be proportional to derivatives of as given by (29), which in general will 
exhibit zero- crossings. 



4 Extensions to the bcisic algorithm 



In this section we briefly describe two extensions to the double null algorithm 
introduced in the previous section — first, to allow for a coordinate system with 
a single null coordinate, and second, extensions to higher dimensional systems. 
We also suggest a technique that can be used to parallelize a characteristic 
code. 

We restrict the discussion on modification for a single null algorithm to the 
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{x,v) = il,0) 



Fig. 6. The fundamental unit cell of the {x, v) coordinate system (2) as depicted in 
Fig. 1, and using the discretization scheme discussed in Sec. 2. Points E,B,C,F 
and D hold the "initial data" for a single evolution step that solves for unknowns at 
point A. The same data structure is used to represent the mesh as with the double 
null scheme (Fig. 3). 

numerical scheme and coordinate system presented in Sec. 2, though in gen- 
eral no significant changes would be needed to alter the algorithm to use an 
outgoing instead of ingoing null coordinate, or use difi^erent FD stencils. 

4-1 AMR with a single null coordinate 

The coordinate system introduced in Sec. 2 has a single null coordinate v 
and a spacelike coordinate x (which becomes timelike when there are trapped 
surfaces). This does not alter the double null algorithm "in spirit", though, as 
demonstrated in Fig. 6 and discussed further below, the spacelike coordinate 
does introduce a preferred direction of integration, and also changes the shape 
of the unit cell in a manner that affects the order in which child cells are 
recursively traversed during evolution. We use the same data structure as 
before to represent the hierarchy, though here north — south links follow lines 
of constant x. We also assume that initial data for the evolution is still specified 
on a pair of null surfaces; this is easy to do in the coordinate system (2) if we 
use a; = 1 as the initial surface of integration in x, for x becomes null in the 
limit x — > 1. 
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Notice from Fig. 6 that because x = const, is timclikc, caiisality forces us ^ to 
integrate in x, along curves v = const., before taking integration steps in v. In 
other words, the causal past of the unknown point A at {x, v) = {xa, va), that 
we want to solve for during a single integration step, now includes regions of 
spacetime with x > xa,v < va and x < xa (hence the extension of the unit 
cell to include point E). Thus, initial data specified along x — xo and v — vo 
will not be sufficient to integrate along the sequence of curves x = Xq — Ax, x ~ 
Xq — 2Ax, instead we need to integrate along v = Vo + Av, v = Vo + 2Av, .... 
Furthermore, note that the size of the integration step Av is now subject to 
a CFL stability condition in that Av must be sufficiently small so that point 
E of the numerical stencil is spacelike separated from A. 

The change in the causal structure of the coordinate system also affects the 
order in which child cells are traversed during the recursive phase of the evolu- 
tion algorithm. For the unit cell depicted in Fig. 6, hnes 18 — 24 of the double 
null algorithm listed in Fig. 4 should be modified to the following: 



18: if (TRE(C)>iiiaxiinum_TRE) then 

19: if (C->child has not been evolved) then evolve_unit_cell (C->child) ; 

20: evolve_unit_cell(C->child->w) ; 

21: evolve_unit_cell(C->child->w->w) ; 

22: evolve_unit_cell(C->child->n) ; 

23: evolve_unit_cell(C->child->n->w) ; 

24: end if 



This sequence of child-cell evolution steps ensures that initial data, consisting 
of points B, C, D, E and F, is always available when we integrate to point A, 
at any level within the hierarchy. Effectively, what we are doing is evolving all 
child points that are contained within the coordinate volume of the unit cell 
to the past of A before evolving to the future of A. The test in hne 19 of the 
modified algorithm prevents the corresponding point from being evolved twice 
in the interior of the grid, which otherwise could happen because neighboring 
unit cells overlap in the x direction. The only place where C — > child is evolved 
is on the initial x = boundary of a refined region or the computational 
domain (where — 1). On such a boundary, if X^ < 1, we initialize the set 
of points corresponding to F in Fig. 6 via interpolation from parent cells (as 
is done with the other points B,C,D and E on the refinement boundary). At 
X.,n = 1, we effectively initialize \1/ at F via extrapolation of the initial data 
ona; = ltoa; = l-|- Ax via = there (the evolution equations for (3 and 
g are first order in x, and do not require initial data at F). 



If we want to maintain a relatively simple time-stepping procedure. 



19 



4.1.1 Initial hierarchy construction and problem dependent options 

Here we briefly describe some features of our one dimensional code that are 
of relevance to a general AMR algorithm, including the interpolation and 
dissipation operators we use, though the particular implementation of these 
features may be problem dependent. 

We have decided to use the function to compute truncation error (TE) 
estimates. This function behaves adequately in tracking regions of high TE in 
the situations we have looked at, except for incoming waves from I~ in the 
vicinity of . The reason why ^ .j, (or any function of ^) fails to give a good 
estimator for the TE of the system there is that, with our choice of coordinates 
and variables, incoming (massless) waves from I~ are propagated essentially 
without error in the vicinity of This has not been a problem for us, as 
the initial data we specify on I~ is always well resolved with a reasonably 
sized coarse mesh. Therefore, the initial hierarchy at a; = 1 only consists of 
two levels — the base level {i = 2) and its shadow {£ = 1). On the other initial 
characteristic v — 0,we generate the initial hierarchy by iterating the following 
until the number of levels stops increasing: we evolve all levels forward one 
coarsest step from v = to v = Avi (refining as usual during the evolution), 
then we reset v to 0, and reinitialize the fields using the hierarchy obtained 
at V = Avi. We start the iteration process with the two coarsest levels fully 
refined. 

Some form of numerical dissipation is usually necessary in Cauchy AMR codes, 
for otherwise the interpolation at grid boundaries is often a source of un- 
wanted, high frequency solution components ("noise"). We have also found it 
necessary to add dissipation to the characteristic algorithm ^ . We do so by 
using a Kreiss-Oliger style filter [31], whereby we modify "^^^^ after it has 
been solved for in (14) via: 

where e should be between and 1 for stability; we typically use e = 0.1 in the 
AMR code. This form of dissipation is adequate in some situations, though 
is not always effective in reducing noise when a new fine level is introduced. 
Furthermore, even though increasing e beyond ~ 0.2 does help to further re- 
duce the high-frequency noise, we have found that it often introduces a small, 
low-frequency, ingoing "tail" component to ^ when the scalar field is predomi- 
nantly outgoing. Thus, work still needs to be done to find a more effective form 
of dissipation for characteristic AMR. A final comment regarding dissipation: 

^ However here the situation is somewhat different, in that we need to apply dissi- 
pation along a "time" direction; in a typical Cauchy AMR codes one only dissipates 
in space. 
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subroutine interpolate (grid function f, point C) 
if (C->parent exists) then 

f [C] = f [C->parent] ; (straight copy) 

else if ( [C->s,P2=C->s->parent,Pl=Pl->s, 
P3=P2->n,P4=P3->n] exist or 
[C->e , P2=C->e->parent , 

Pl=Pl->e,P3=P2->w,P4=P3->w] exist ) then 

f[C] = 1/16* (-f [PI] + 9*(f [P2]+f [P3]) - f[P4]); 
('centered' cubic polynomial interpolation) 

else if ( [C->s,P3=C->s->parent,P2=P3->s, 
Pl=P2->s,P4=P3->n] exist or 
[C->e , P3=C->e->parent , P2=P3->e , 
Pl=P2->e,P4=P3->w] exist ) then 

f[C] = 1/16* (f [PI] + 5*(f [P4]-f [P2]) + 15*f[P3]); 
('backwards' cubic polynomial interpolation) 

else if ( [C->s,C->s->e,PA=C->s->e->parent, 

PB=PA->n,PC=PA->w,PD=PC->n] exist or 
[C->w , C->w->s , PA=C->w->s->parent , 
PB=PA->n,PC=PA->e,PD=PC->n] exist ) then 

f[C] = 1/4* (f [PA] + f[PB] + f[PC] + f [PD]); 
(bilinear interpolation) 

else if ( [C->s,Pl=C->s->parent,P2=Pl->n] exist or 
[C->e,Pl=C->e->parent,P2=Pl->w] exist or 
[C->w,Pl=C->w->parent,P2=Pl->e] exist ) then 

f [C]=l/2*(f [PI] + f[P2]); (linear interpolation) 

end if 

end of subroutine interpolate 

Fig. 7. A pseudo-code description of the interpolation routine we use to initialize 
functions at points on the initial data surfaces of interior fine levels. We use cubic 
interpolation if it is possible to do so given the structure of adjacent points on 
the parent level, otherwise we switch to lower order interpolation. For the kinds 
of hierarchies generated by the algorithm described here, one of these conditions 
will always be matched, hence the "if" part of the last "else if" statement is not 
necessary. 

after a new fine level is introduced, we typically disable regridding in a small 
buffer zone next to the boundary of the new level; this gives the dissipation 
some time to work at eliminating the noise introduced via interpolation, and 
prevents a refinement "cascade" . 

We use cubic Lagrange interpolating polynomials to initialize fine level grid 
functions at refinement boundaries if a sufficient number of adjacent points 

are available on the parent level, otherwise we use linear interpolation — Fig. 7 
is a pseudo-code description of our interpolation routine. We have also exper- 
imented with linear interpolation for all child points, though found that this 
produces significantly more noise after refinement. 
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4-2 Application to higher dimensional systems 

The extension of the AMR algorithm to higher dimensions, in other words 
to problems where there is dependence on additional spacelike coordinates z\ 
i — 1,2, ..,d, is straight-forward. The coordinates are treated just as the 
spacelike coordinates are in the traditional BhO algorithm. The null AMR 
evolution algorithms described in the preceding sections are un-modified, ex- 
cept that now at each "point" of the mesh structure one needs to store a 
list of d-dimensional arrays. Thus, the composite grid hierarchy at any point 
{u,v) (or {x,v), etc.) within the computational domain will look like a B&O 
hierarchy for a (i-dimensional problem. The particular set of arrays needed at 
a given point and at a given level in the hierarchy are determined, as usual, by 
local truncation error estimates. Consequently, standard clustering algorithms 
will be needed to convert the region of high TE to a set of grids. 



4.-3 Parallelization 

Here we briefly mention a scheme that could be used to parallelize a charac- 
teristic code, with AMR or otherwise. For simplicity, we illustrate the concept 
with a unigrid double null scheme, though it is not difficult to generalize it. 
The idea is to use a set of n processors as a pipeline, as demonstrated in 
Fig. 8 below. Integration starts on a single node, where processor 1 solves the 
equations within a region Ri of size A — AU AV to the future of the initial 
point {uq.Vq). Afterward, initial data is available to simultaneously solve the 
equations on two regions R2 and R-^, of the same size A, to the future of Ri. 
Processor 1 (arbitrarily) proceeds to solve the equations in region R2, while 
supplying the relevant initial data to processor 2 to evolve region R3. After 
processors 1 and 2 have solved the equations in regions R2 and R3 respectively, 
initial data is now available to simultaneously evolve three regions of size A. 
Processor 3 can now enter the pipeline, and the process continues. After the 
first n steps, the pipeline will be full, and all nodes will be involved in the 
computation. 

One desirable feature of this algorithm (compared to a typical parallelization 
scheme for a Cauchy problem) is that it would be possible to stagger the 
communication — in other words, all nodes do not need to, and ideally should 
not, communicate at the same time. Furthermore, in a certain sense the com- 
munication is one way, so that when a processor is finished solving one block 
of data it can simultaneously send initial data to the next processor down 
the line while starting to evolve a new block. It would also not be difficult to 
modify the algorithm to dynamically subdivide the region of space a processor 
must solve in a single step, to provide better load balancing among the nodes 
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Fig. 8. An illustration of the technique suggested to parallelize a characteristic 
code. In this example, three processors are available to simultaneously solve the 
equations. We subdivide the grid into blocks of equal area, and assign blocks to 
individual processors as shown in the figure. The heavy line drawn on the grid at 
step i represents the surface where initial data will be available at step i + Before 
step one, there is only a single block where sufficient initial data is available to solve 
the equations, hence only one processor is active then. After step one there are two 
unsolved blocks that can be evolved, and thus two processors become active; and 
so on, until the "pipeline" of processors is full. 

(this will be necessary when AMR is used, or in a higher dimensional simula- 
tion if the total number of points on grids in the extra spacelike dimensions 
depend on where in {u,v) space one is). 



5 Tests & Results 



In this section we present results from two evolutions obtained using the ID 
characteristic AMR code described in the previous sections. The first example 
models a black hole interacting with an initially outgoing pulse of massive 
scalar field energy. The relatively large mass term we use, in conjunction with 
the compactified radial x coordinate, causes a wide range of relevant scales to 
develop at late times, and we show that the algorithm is able to track these 
features with reasonable accuracy via comparison with results from unigrid 
evolution. The second example shows a near-critical collapse [32] of the mass- 
less scalar field. For reasons we will explain below, our coordinate system is 
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not well suited to studying this kind of critical phenomena, and hence the 

example is not very close to criticality. Nevertheless, we are close enough to 
demonstrate that the algorithm can also adapt to the very small length scales 
that develop from ingoing initial data. 



5. 1 Massive scalar field - black hole interaction 



For the first example, we specify ^ along i> = as follows: 



'^{v — Q,x)— Ax{l — — xjx-if ^ X\ < X < X2i 

— elsewhere, (31) 

choose Xx — 0.15, X2 = 0.25, A — 5x10^ and the mass parameter m — 5. We 
set g{x — l,v — 0) — 0.2, so that the asymptotic mass of the spacetime is 0.1; 
the initial {v = 0) black hole mass is ~ 0.063, indicating that the scalar field 
contributes ~ 0.037 to the initial mass of the spacetime. We ran the simulation 
to f = 10, by which time the black hole mass had grown to roughly 0.068 due 
to accretion of scalar field energy. First, to demonstrate the correctness of our 
solution to the finite difference equations, in Fig. 9 we show the convergence 
factor for the scalar field 

11*4/1 - *2/i lb /oo^ 

\\^2h — ^h\\2 



computed using solutions from 3 different resolutions of unigrid simulations. 
*4ft has 1025 points in x, *2/i 2049, and 4097. The Courant factor is 0.25 
in all cases; i.e. Av = 0.25Aa;. That is around 4 for most of the simulation 
indicates that we are seeing second order convergence, as expected (and that 
Q-q, eventually starts to deviate from 4 is also not unexpected — accumulating 
numerical errors, especially in our compactified coordinate system, will even- 
tually cause any finite resolution simulation to move away from the convergent 
regime). For brevity we do not show convergence factors for the other fields; 
they also exhibit second order convergence. To demonstrate the accuracy of 
the adaptive code, we first show a comparson between a solution generated 
with AMR to the finest unigrid solution, and then present some results from 
a convergence test. 



5.1.1 Comparison between AMR and unigrid solutions 

For the comparison between the AMR and unigrid results, the parameters 
for the AMR run were set so that the base level (2) has 513 points (hence 
the shadow level (1) has 257 points), and the maximum allowed TE is such 
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that early on (in v) during the evolution the finest level is 5, giving the same 
resolution locally as that of the finest unigrid simulation. At late times during 
the adaptive run additional levels are introduced to track the outgoing pulse 
(whose width shrinks due to the radial coordinate we use) — see Fig. 10 for 
several snapshots of ^ along v — const, surfaces during evolution, Fig. 11 for 
the maximum level as a function of x at the same instants of v shown in Fig. 
10, and Fig. 12 for a plot of the maximum level in the hierarchy over all a; as a 
function of v. To gauge the quality of the adaptive solution, we use the highest 
resolution unigrid solution for ^ as a benchmark, and compare this solution 
to the lower resolution unigrid and adaptive results. Fig. 13 shows the £2 norm 
(computed along v = const, slices of the spacetime) of these differences, and 
Fig. 14 plots \t'(x, f = 10) from the four simulations near the leading edge 
of the pulse, for visual comparison. Fig. 13 demonstrates that early on the 
adaptive solution gives slightly worse results than the lower resolution unigrid 
solution, which is not too surprising as the AMR solution only covers a portion 
of the computational domain with comparable or higher resolution. However, 
at late times the adaptive solution starts to outperform the coarser unigrid 
solutions as the AMR is able to keep the narrowing pulse well resolved. 

5.1.2 AMR convergence test 

For a separate convergence test of the AMR code we ran three simulations 
with the same initial data as described above, varying the base grid resolution 
and maximum allowed truncation error to mimic doubling the resolution from 
one run to the next. Specifically, the lower resolution simulation had 257 points 
in the base level (2) with a maximum TE of r^h, the medium resolution run 
had a base grid of 513 points and a maximum TE of r2h = T4/j/4, while the 
higher resolution simulation had a base grid of 1025 points and a maximum 
TE of Th — T4h/i^- However, due to limited available computer resources we 
restricted the maximum depth of the hierarchy to 7 for these three runs (at this 
stage our code is not memory efficient, and the high resolution run was using 
most of the available memory nearing the end of the simulation) . Note that this 
scheme for convergence testing an AMR code will not produce identical grid 
hierarchies that differ only in resolution, because the truncation error estimate 
used in each numerical simulation will not scale exactly as the leading order 
part of the actual truncation error, which decreases by a factor of 4 each time 
the mesh spacing is halved for a second order accurate scheme. Nevertheless, 
if the test were to show non- convergent results it would be a clear indication of 
problems with the implementation. Fig. 15 below shows the convergence factor 
for ^ (where to calculate Qq, as in (32) we first interpolated the solutions to 
identical uniform grids). Early on we get a convergence factor that is closer to 
first than second order convergence. The primary reason for this appears to be 
grid-boundary "noise" generated at parent-child interfaces, and as discussed in 
Sec. 4.1.1 our current interpolation/dissipation scheme is not yet very effective 
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at reducing this noise. Later on in the simulation the convergence behavior 
appears to improve significantly (and become unrealistically high), however 
this is mostly because we had to limit the maximum level to 7. At late times 
this is not sufficient to maintain the desired TE (see Fig. 11 for the level 
structure generated by the AMR run described in the previous section, which 
had a base resolution equivalent to that of the medium resolution simulation 
here), and the solutions start to drift away from the convergent regime. In fact, 
by V = 10 the lower resolution run had accumulated significant phase errors 
in whereas the medium and high resolution solutions were still roughly in 
phase, which gives some explanation for the anomalously high value of 
(32) in this case. 



5.2 Massless scalar field critical collapse 

For a second, brief example, we consider the near-critical collapse of an initially 
ingoing pulse of the massless scalar field: 



and set Vi — 0.01, V2 — 0.11 and g{x — l.,v — — A \s tuned (via a 
bisection search) so that the collapsing pulse is close to the threshold of black 
hole formation. As mentioned before, the coordinates we use are not well 
suited to studying type 11 critical collapse, where one should be able to form 
arbitrarily small black holes in the super-critical regime. The reason is that in 
our coordinate system, a non-zero initial value g{x = 1,^ = 0) = 2mo describes 
a spacetime containing a black hole of mass mo (assuming ^(x, v = 0) = 0). 
Truncation error effects in the integration oi g{x = l,v), and subsequent 
evolution in x, causes small errors in g that effectively behave as if a small 
black hole (of size proportional to the TE) had been present in the initial 
data. This is not a problem for unigrid evolution, as the erroneous black hole 
is typically smaller than the grid spacing; however during a critical evolution 
where arbitrarily small scales unfold, and refinement resolves these scales, this 
black hole is eventually revealed. Thus the resolution of the initial data at /~ 
places a limit on how close to critical we can evolve ^ . For this example we 
chose a base (level 2) resolution of 2049 in x (with Av — 0.25 Ax); then the 
smallest black hole we can form is on the order of 10~^, which is roughly 1/2 
the size of a base-level cell in x. Fig. 16 shows several snapshots of ^/r from 



To seriously study critical collapse with this characteristic AMR algorithm one 
would therefore need to choose more appropriate coordinates, such as one based on 
an outgoing null coordinate for instance [33]. 



1)^A{1-V/Vi)'^{l-V/V2)^, Vi<V<V2, 



— elsewhere. 



(33) 



26 



1 — I — r 




J \ L 



2 4 6 8 10 

V 



Fig. 9. The convergence factor (32) from unigrid evolution of the black 
hole-massive scalar field initial data. 

the evolution of the nearest-to-critical sub-critical amplitude we found, and 
Fig. 17 shows the corresponding level structure. 



6 Conclusion 

In this paper we have introduced a new algorithm that can be used to add an 

adaptive mesh refinement framework to a characteristic evolution code. The 
algorithm is similar to the Berger and Oliger algorithm for Cauchy codes, and 
shares many of its desirable features, including dynamical regridding via local 
truncation error estimates, efficient use of computational resources via the re- 
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Fig. 10. ^{x) along several v = const, slices of the spacctimc, from the adaptive 
black hole-massive scalar field simulation discussed in section 5.1.1. The position of 
the excision surface is denoted by a vertical dashed line, and is always kept a distance 
of 0.025 in x inside the apparent horizon; the scalar field is set to zero inside the 
excision surface. The scalar field pulse is initially outgoing, and the higher-frequency 
components of the field continue to travel outward at the speed of light. The non-zero 
mass term in the scalar field equation of motion causes lower frequency components 
of the field to travel at speeds less that 1; these trail the leading pulse of the field, 
and interact more strongly with the black hole. See Fig. 11 below for a plot of the 
maximum AMR level along the same slices shown here to see how the AMR algo- 
rithm correctly tracks the outgoing, higher-frequency components of the solution. 
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Fig. 11. Maximum level of the hierarchy from the adaptive black hole-massive scalar 
field simulation discussed in section 5.1.1, at the same slices of the computational 
domain shown in Fig. 10. 



cursive evolution scheme, and in principle it does not require modifications to 
the underlying unigrid finite difference scheme. As discussed in the introduc- 
tion, we believe AMR is essential to achieve accurate results from simulations 
in many general relativistic scenarios. Based upon the early success of this 
AMR technique with the one dimensional code presented here, we think it 
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Fig. 12. Maximum level of the adaptive hierarchy along v = const, slices of the 
spacetime, from the adaptive black hole-massive scalar field simulation discussed in 
section 5.1.1. 

would be well worth the effort to apply the method to higher dimensional 
problems of interest. 
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